Critical ionic transport across an oxygen-vacancy ordering transition

Phase transition points can be used to critically reduce the ionic migration activation energy, which is important for realizing high-performance electrolytes at low temperatures. Here, we demonstrate a route toward low-temperature thermionic conduction in solids, by exploiting the critically lowered activation energy associated with oxygen transport in Ca-substituted bismuth ferrite (Bi1-xCaxFeO3-δ) films. Our demonstration relies on the finding that a compositional phase transition occurs by varying Ca doping ratio across xCa ≃ 0.45 between two structural phases with oxygen-vacancy channel ordering along <100> or <110> crystal axis, respectively. Regardless of the atomic-scale irregularity in defect distribution at the doping ratio, the activation energy is largely suppressed to 0.43 eV, compared with ~0.9 eV measured in otherwise rigid phases. From first-principles calculations, we propose that the effective short-range attraction between two positively charged oxygen vacancies sharing lattice deformation not only forms the defect orders but also suppresses the activation energy through concerted hopping.


Supplementary Figures
Supplementary Fig. 1: Structures of as-grown BCFO films on SrTiO3 substrates. a, The surface topographic images of 100-nm-thick as-grown BCFO films on SrTiO3 substrates have atomically flat surfaces. The surface roughnesses were evaluated to be less than 0.3 nm. The film of xCa = 0.1 has a monoclinic structure with ferroelectricity similar to BiFeO3. Other BCFO films at xCa = 0.2 or higher have pseudo-tetragonal structures. Scale bars represent 1 m. b, X-ray 2θ-ω scans at room T for samples with varying xCa (0.1 ~ 0.6). Dashed line represents the (001) peak of SrTiO3. BCFO films at xCa = 0.2 or above have superlattice peaks (green lines, SL±), indicating the periodic appearance of VO layers along the c axis. c, c-axis lattice parameter versus xCa. The red line denotes the linear interpolation between the lattice parameters of BiFeO3 and CaFeO2.5 according to Vegard's law. Error bars were estimated to be 1/5 of the FWHMs of the BCFO film peaks. d, Periodicity of ordered VO layers versus xCa. Error bars were estimated from the FWHMs of the SL± peaks. Red line denotes an empirical rule, i.e., nperiod = 1.5/xCa. As the distance between VO layers becomes closer, it tends to deviate from the empirical rule. BCFO films of xCa = 0.2 (0.3) have nperiod = 7 or 8 (5 or 6) unit cells, but for xCa = 0.6, most of the VO layers are ordered with a periodicity of 3 unit cells. where Re(ion) is electronic (ionic) resistance. In the low-frequency regime, the total resistance (Rbulk + Rt) can be interpreted as Re due to ion-blocking electrodes. [100] [010] [100] [001] c d e f a especially out-of-plane oxygens, are found to be the most prominent around the clustered VOs. This indicates that in the clustered case the relaxation energy gain was achieved mainly by oxygen ions and was large enough to overcome the repulsive Coulomb interaction between VO 2+ s. c-f, Atomic structures of four period types of VO layer in a representative BCFO (xCa = 0.25). c, 2-, d, 4-, e, 6-, and f, 8-unit cell periods with the same xCa and δ are, respectively, constructed for nperiod = 0.5/xCa, 1.0/xCa, 1.5/xCa, and 2.0/xCa with 1, 2, 3, and 4 VOs in the VO layer.

Supplementary Note 1: A plot of position-time curves on a linear time scale
We provide the same plot of the position-time curves in Fig. 1c on a linear scale for comparison ( Supplementary Fig. 11). The plot on a linear time scale shows the low-temperature region clearly, while the plot on a logarithmic time scale displays the high-temperature region more clearly. The experimental data match well with the model fit function except for the initial time regime. The measured phase boundary near the electrode tends to be more slowly moving than the theoretical curve. The nucleation of dark conducting phase takes time after voltage is turned on and this time is now considered as a fitting parameter of t0 in the model. But the physical origins of determining t0 are hardly understood yet. There are uncharacterized effects such as the interfacial effect between Pt metal and BCFO film, the extra BCFO region underneath the metal electrode, and the nucleation times for intermediate and dark phase evolution. The fitting was made for a well-defined region of the dark-phase propagation between ~50 µm and 320 µm where the interfacial effect near the electrode is thought to be neglected.

Supplementary Note 2: A more general form of the Nernst-Einstein relation
A more general form of the Nernst-Einstein equation is discussed to describe the diffusion phenomena in non-dilute systems. Before going into details, we need to clarify two different concepts of diffusion coefficients. The chemical diffusion coefficient (̃) is a proportionality constant between the gradient in concentration and the diffusion flux described from the Fick's law. On the other hand, the conduction diffusion coefficient ( ) is defined from the classical Nernst-Einstein equation ( = B / , with conductivity , Boltzmann constant kB, temperature T, concentration n and charge q). From the Green-Kubo relation, they are defined as [Gomer, R. Rep. Prog. Phys. 53, 917 (1990)]: N corresponds to the total number of particles, ( ) is the fluctuation in the particle number, its ensemble average value has a relation of 〈( ) 〉 = 〈 〉 − 〈 〉 . Δ is the displacement of i th ion in d-dimensional system after time t, and the time-normalized vector corresponds to the velocity of i th ion. The conventional approaches such as the electrochemical impedance spectroscopy, galvanostatic (or potentiostatic) intermittent titration techniques, and even the electrocoloration method measure the chemical diffusion coefficient ̃ in terms of the fact that a potential gradient is introduced.
On the other hand, the tracer diffusion coefficient ( * ) is defined as It carries the elementary picture of diffusion for a single particle executing a random walk in ddimensional space. We notice the correlation terms between different particles are omitted, in contrast to the ̃ and . Isotope diffusion experiments using less concentrated ad-particles can give the value of * .
The * can be related to by multiplying the correlation factor : (4) Provided that a particle is weakly correlated with other particles such as in dilute systems, i.e., →∞ ⟨ ( ) • ( )⟩ → ( ≠ ), * becomes equal to , i.e., → 1. Observing the deviation of from 1 by performing an isotope diffusion experiment would be an interesting experiment to test the correlation effect.
Since our experiments are the electrochemical impedance spectroscopy and the electrocoloration approach, we need to focus on the relation between ̃ and . As shown in the work by G.  1569-1578 (1977)]. Of course, it also approaches to 1 in dilute systems. In the as-grown samples without nitrogen annealing, the BCFO samples are mixed ionic-electronic conductors. Also, the electronic conductivity predominates the ionic conductivity because the electronic carriers have a much greater mobility than the ionic species. Since electronic carriers are generated with oxygen vacancies, the electroneutrality condition i ∆ i ≈ ∆ e is satisfied in our system ( i is the charge number of an ionic defect, i and e are the concentrations of ionic and electronic carriers). In the constraints together with obeying the Raoult's law, the enhancement factor can be written as, The electrocoloration method quantifies the mobility and thus the diffusivity deduced by the Nernst-Einstein relation corresponds to . On the other hand, the electrochemical impedance spectroscopy is known as a technique for characterizing ̃ in the electrolyte regime. By comparing the results of the two methods, we can determine the value of to be ~2 (See Fig. 1f and Supplementary Fig. 12c). Considering a VO donates two electrons, e is estimated to be a quarter of i . Moreover, is just a weighting factor independent of T in the experimental range and thus the overall tendency for the EA is not changed.
As abovementioned, the exact value of can be empirically determined on the assumption that ionic defect density of i is appropriately specified, but the estimate is not simple because the spatial distribution of VOs in BCFO is not fully randomized but they partially play a role of constructing an ordered structure of VO channels. Taking into account a deviation () of oxygen stoichiometry from the ideally compensated semiconductor, Bi1-xCaxFeO3-x/2+, we may regard the concentration of only the excess oxygen ions as i . Since a single excess oxygen ion is an acceptor for generating two holes ( i = ) and thus the concentration of produced electronic hole carriers is e = | i | i , the value of becomes 3.
It is worthwhile mentioning that the conductivity of color tracking method could be overestimated by a factor of 3/2. We determined the ionic conductivity = from the measured mobility. When the propagation of the phase boundary between the dark conducting phase (phase III with 0 VO) and the dark-yellow phase (phase II with 2 VOs) is examined, the number of moving VOs are only 2 VOs per 2×2 pseudocubic unit cells in an oxygen-deficient layer. As compared with the as-grown phase (phase I with 3 VOs), the change in ionic carrier density is two thirds of and so does . If we accept for the ~3 , the conductivities from the electrocoloration method and the impedance spectroscopy ( Supplementary Fig. 12) are exactly matched each other relying on the correction.

Supplementary Note 3: Nitrogen annealing effect
A previous paper [N. Masó and A. R. West, Chem. Mater. 24, 2127-2132 (2012)] motivates us to anneal our samples in the N2 gas environment. The paper has reported that the exact stoichiometric Bi0.7Ca0.3FeO2.85 ceramics exhibit purely ionic conduction based on impedance spectroscopy. Small amount of off-stoichiometric defects/ions have largely changed conductivity and activation energy barrier. In the Masó and West paper, annealing under oxygen or air environment slightly changed the weight of crystal relative to the nitrogen annealed one, i.e., oxygen annealing: Bi0.7Ca0.3FeO2.85+0.016, and air annealing: Bi0.7Ca0.3FeO2.85+0.01. Even if the variations look tiny, those give rise to large differences in conductivity and activation energy, indicating an electronic-carrier dominant conduction.
To compare the video method and impedance measurement, we need to make ionic-conduction dominant Bi0.7Ca0.3FeO2.85 thin films to determine ionic conductivity by using the standard AC impedance measurement. The as-grown thin film shows an electronic-conduction dominant state without any Warburg-spike feature (red data in Supplementary Fig. 12a). We put a lot of effort into finding the optimal nitrogen-gas annealing conditions for making the ion-conduction dominant state in our Bi0.7Ca0.3FeO2.85 films. We successfully obtained the Bi0.7Ca0.3FeO2.85 films by annealing at N2 gas pressure of 10 mTorr and 650 °C for 24 hours to suppress electronic carriers induced by offstoichiometric oxygen atoms. We used platinum interdigitated electrodes to apply an AC bias of 100 mV with a frequency sweep ranging from 1 MHz to 0.1 Hz. We observed a large semicircle and a lowfrequency branch (black data in Supplementary Fig. 12a). The large semicircle is attributed from the ionic conduction in bulk, and the branch in the extended low frequency part is a feature of the wellknown Warburg spike, arising from the existence of a constant phase element often found in ionic conductors [J. Jamnik and J. Maier, Treatment of the impedance of mixed conductors equivalent circuit model and explicit approximate solutions, J. Electrochem. Soc. 146, 4183-4188 (1999)].
The diameter of the large semicircle represents the ionic resistance of Bi0.7Ca0.3FeO2.85. Impedance measurements at different temperatures also show a similar Warburg-diffusion behavior and enable us to estimate the temperature-dependent ionic conductivity ( Supplementary Fig. 12b). We evaluated ionic conductivity of Bi0.7Ca0.3FeO2.85 thin films using the video analysis by tracking the motion of phase boundary [J. S. Lim et al., Ultrafast collective oxygen-vacancy flow in Ca-doped BiFeO3, NPG Asia Mater. 2, 084412 (2018)]. Comparing the ionic conductivities obtained by two different methods, they match well with each other, revealing the same order of magnitude of conductivity and the same thermal activation energy (Supplementary Fig. 12c). The sample used for the video recording method was not annealed. Nonetheless, the mutual consistency suggests that the visualization technique provides a useful pathway for characterizing ionic conduction properties regardless of the presence of significant electronic conduction. Supplementary Fig. 12: Impedance measurements of Bi0.7Ca0.3FeO2.85+. a, Nyquist plots of asgrown and N2 gas annealed Bi0.7Ca0.3FeO2.85+ at 300 ºC. b, Nyquist plots of N2 gas annealed Bi0.7Ca0.3FeO2.85 with respect to temperature. c, A comparison of ionic conductivities obtained from the two different methods, video analysis (orange balls) and impedance measurement (black balls).

Supplementary Note 4: Determination of the oxygen tracer diffusivity
Tracer diffusivity (D*) of oxygen ions in a BCFO thin film (xCa = 0.45) was investigated by an oxygen isotope exchange profiling technique using time-of-flight secondary ion mass spectroscopy (TOF-SIMS). A ~50 nm-thick LaAlO3 capping layer was deposited on the ~55 nm-thick BCFO thin film at xCa=0.45 to determine the tracer diffusivity along the in-plane direction. The sample was scratched with a straight line by a scalpel to open a lateral edge of the thin film to an atmosphere ( Supplementary Fig.  4a). Thus, the oxygen exchange was allowed along the in-plane direction [100], while prevented along the out-of-plane direction by the capping layer. The topographic image in Supplementary Fig. 4b shows that the gap was clearly made along the [010]. The width and depth of the gap was ~7 μm and ~200 nm, respectively. The sample was evacuated to <5×10 -3 Torr and then pre-annealed at 430 °C for 1 hour under a stream of pure oxygen gas (99.995%) of natural isotopic abundance. After the pre-annealing process, the sample was cooled to room temperature and evacuated again to <5×10 -3 Torr. Then, 18 O2rich gas (97.15%) was inserted up to a pressure of ~500 Torr, and annealing at 350℃ in the isotope gas environment was conducted for 90 minutes. At last, the sample was cooled down to room temperature.
I( 18 O -) is the 18 Ointensity, IOis the sum of the 18 Oand 16 Ointensities, Cbg is the natural isotopic fraction of 18 O, Cg is the isotopic fraction of the gas used for the annealing process, D* is the bulk oxygen tracer diffusivity, and k is the tracer surface exchange coefficient. Cbg and Cg were set to 0.002 and 1. In our analyses, IOwas also set as a fitting parameter. The data points near the scratch (closer than 0.003 cm) were excluded for the fitting to avoid any unwanted edge effects and to focus on pure bulk oxygen diffusion. The deviation of the FeOintensities at the edge region also suggests the existence of an edge effect.
As a result, D* = 2.77(± 0.19)×10 -9 cm 2 /s and k = 1.46(± 0.04)×10 -9 cm/s were determined at 350 °C in the BCFO film of xCa = 0.45. The determined diffusivity value is approximately twice lower than the value of conduction diffusivity ( ) of 5.9×10 -9 cm 2 /s at the same composition and temperature. A slightly underestimated value seems reasonable as tracer diffusivity does not account for correlation effects due to dilute concentrations of isotope ions.
Depth profiles of 18 Ointensity through the entire LaAlO3 capping layer and the upper BCFO part were analyzed ( Supplementary Fig. 4e) in three regions with different distances from the gap. The red points represent the average depth profile in the area R1 (shown in the red rectangle in Supplementary Fig. 4c) 30-50 μm distant from the edge of gap. Similarly, the green (blue) points are average depth profiles in the area R2 (R3) 50-70 (70-90) μm distant from the edge of gap. The 18 Ointensity abruptly decays within a sputtering time of 200 s. Considering that ~2000 s was taken to fully etch the 50 nm-thick LaAlO3 layer, the sputtering time of 200 s corresponds to ~5 nm. Thus, most of 18 O ions injected to LaAlO3 could not diffuse deeper than 5 nm in all the R1, R2, and R3 regions. By etching longer than 2000 s, we can remove the 50-nm-thick capping layer completely and reach the BCFO thin film. In contrast to the almost same profiles of the LaAlO3 capping layer, a gradual decrease in 18 Ointensity can be observed from R1 to R2 to R3 regions in the BCFO layer, indicating that the concentration gradient was created by oxygen tracer exchange and diffusion through the in-plane direction.

Supplementary Note 6: Theoretical understanding of defect orderings and competition
The atomic structure of as-grown VO-ordered BCFO superlattice is modeled based on two known conditions: (i) the VO concentration (δ) should follow δ = xCa/2 because VO in BFO is a double-electron donor (VO 2+ ) compensated by two CaBi 1acceptors; (ii) the distance between neighboring oxygendeficient layers (nperiod) is experimentally measured to be nperiod = 1.5/xCa. To meet these conditions, three VO 2+ s should be present per 2×2 cubic perovskite unit cell with eight oxygen sites (see Supplementary  Fig. 9a for a representative BCFO atomic structure with xCa = 0.25 and nperiod = 6, based on the calculational results below).
To understand the microscopic origins of superlattice formation in BCFO films, we performed firstprinciples DFT calculations. Regarding the mysteriously ordered VO layers in BCFO, we identified the following three questions for theoretical investigation: (1) Why does VO prefer to cluster together and form the ordered layers? (2) Why does the nperiod follow the rule nperiod = 1.5/xCa? (3) What are the atomic models for the [100] and [110] patterns of the VO channels? By answering these questions, we might be able to better understand the experimental findings.
(1) Ordered VO cluster formation: We first studied various configurations of two VO 2+ defects. In our model calculations (see Methods of calculation), the VO 2+ -pair configurations at and near the closest distance are found to be more stable than the others. This pairing trend is not common in terms of Coulomb interaction between defects. Thus, we examined the origin of pairing through the atomic and electronic structures. We found that, without atomic relaxation, VO 2+ defects prefer to be separated ( Supplementary Fig. 9b), in agreement with their repulsive Coulomb interaction. When we allowed atomic position relaxation, however, the 2VO 2+ cluster was favored over the separated one by 0.86 eV. This means that the atomic relaxation energy in the 2VO 2+ cluster is large enough to overcome the repulsive Coulomb interaction between VO 2+ s ( Supplementary Fig. 9b). Our detailed analysis shows that the out-of-plane apical oxygens in the 2VO 2+ cluster undergo the most pronounced atomic relaxation (Supplementary Fig. 10a and b). This can be understood with crystal field splitting of Fe orbitals in the polyhedra near the 2VO 2+ . While a pair of VO 2+ form two five-fold coordinated square pyramids, the 2VO 2+ cluster naturally relaxes to a four-fold coordinated tetrahedron adjacent to two square pyramids. With the crystal field splitting of the Fe-d orbitals, the relaxation of the tetrahedron requires atomic displacements large enough to reverse the order of t2 and e bands, whereas the relaxation of the square pyramid only occurs to a small extent with no significant reordering of Fe d orbitals (t2g  dxy and dxz/dyz, and eg  d 3z 2 −r 2 and d x 2 −y 2). Our electronic structure analysis demonstrates that the Fe ion centered at the tetrahedron obtains a large energy gain through the atomic relaxation, in contrast to the square pyramid ( Supplementary Fig. 9c). Moreover, the square plane is turned out to be highly unstable and thus we consider only three classes of 3VO configurations that are composed of two tetrahedra and two square pyramids, as shown in Supplementary Fig. 9d. These were evaluated to be stable (at least 42-154 meV/f.u.) compared to the other configurations.
(2) nperiod = 1.5/xCa rule: To evaluate the most favorable value of nperiod in a representative BCFO at xCa = 0.25 and δ = xCa/2, we considered 2-, 4-, 6-, and 8-unit-cell periods, corresponding to nperiod = 0.5/xCa, 1.0/xCa, 1.5/xCa, and 2.0/xCa, respectively (Supplementary Fig. 10c-f). Since the total number of VOs is conserved, the VO density within the deficient layer increases accordingly. Supplementary Fig. 9e shows that the 6-unit-cell period corresponding to nperiod = 1.5/xCa, which removes 3 out of the 8 oxygens in the oxygen-deficient layer, is the most stable one as a result of the optimization of atomic relaxation energy gain (bottom of Supplementary Fig. 9e).
According to HK maps in Supplementary Fig. 8a, all films are fully strained and the in-plane lattice parameters are identical to those of substrate. The average c-axis lattice parameter (cavg) was obtained from the lattice parameter of the supercell formed by the VO ordering. We also considered local c-axis lattice expansion at oxygen-deficient layers, while maintaining the cavg value (smaller than the expanded c) for the other layers. The c-axis lattice parameters were decided based on local tetragonal c/a distortion values ( Supplementary Fig. 5) and by comparing with the experimental X-ray diffraction data ( Supplementary Fig. 8). In the L scans at half-ordered in-plane positions, additional peaks appear in the middle of BCFO film peaks in (0 0 L) scans, which could be explained by domains with VO channels stacked in an out-of-phase fashion with an in-plane one-unit-cell shift. The relative populations of VO domains were obtained by the comparison between the structure factor calculation with spatial coherence of possible domain distributions and the intensity modulation in L scans. The resultant parameters and possible domains are summarized in the Supplementary Table 1. A spatial coherence between multiple VO domains arising from different stacking orders of the oxygen deficient layers is essential for explaining the intensity modulation in the L scans ( Supplementary Fig. 15).
Supplementary Fig. 15: Spatial coherence between domains with different stacking order of the oxygen deficient layers. a, Monotonic change appears in the case of spatial incoherence. b, Relative peak intensity become altered in the case of spatial coherence. Here, I, ai and Fi stand for intensity of diffraction, population of ith domain, and structure factor of ith domain. The red dotted circle and red dotted lines indicate VO and oxygen deficient layers.  Table 1: Parameter values for structure factor calculations. The parameter values and the relative domain population used in structure factor calculations. The cVo, co and cavr correspond to c-axis lattice parameters of oxygen-deficient layers, the layers without VO channels and their average. Each relative population of domains corresponds to the inset unit-cell schematics above each value. The red dotted lines indicate oxygen-deficient layers and the blue dotted lines indicate the ones with the inplane one-unit-cell shift.

Supplementary Note 8: Analysis of the propagation of intermediate phase region
We have observed three color phases in the color diagram: pale yellow (as-grown), dark yellow (intermediate) and dark black (completely formed) phases. We have mainly dealt with the distinct color boundary between dark-black and dark-yellow phases to characterize the motion of oxygen vacancies in the main manuscript. Here, we analyze the front propagation of the intermediate phase as well. Supplementary Fig. 16: Temperature dependence of the phase front mobility for the intermediate and dark phase propagations. a, Arrhenius plots of the product of collective VO mobility (µ) and temperature (T) versus reciprocal temperature for different xCa = 0.2 ~ 0.5. b, The activation barrier EA with respect to xCa, which is obtained by the slopes of linear fit lines in a. Solid and open symbols (solid and dashed lines) correspond to the intermediate phase and the dark conducting phase, respectively. We note the boundary between the intermediate phase and the pristine phase is not clearly seen in BCFO films at xCa = 0.1 or 0.6.
The color boundaries between the intermediate dark-yellow phase and the original yellow phase at different temperatures are highlighted by green dashed lines in the Supplementary Fig. 2. Each of the lines was obtained by fitting the position-time (xint vs. tint) curve of the boundary of intermediate phase to the same equation which was used for the analysis of dark-phase propagation. We can estimate the ionic mobility µ from the experimentally determined Δtint=tint -t0 through the relation µ = L 2 /(2ΔtintV). L stands for the channel length of 400 µm and V is the applied voltage. The error bar of Δtint was set to be 20 % of the corresponding value. The ionic mobility of intermediate phase propagation shows a similar value compared to that of the dark-phase propagation at the examined xCas ( Supplementary Fig.  16a). The product of temperature and mobility is also linearly proportional to inverse of temperature, indicating EA from the slope. The EA of ionic transport for the intermediate phase is roughly ~ 0.1 eV higher than that of the dark phase ( Supplementary Fig. 16b). The value of EA can be larger as much as ~0.4 eV in the case of individual hopping, according to the density functional calculation result. Accordingly, it suggests the correlation effect cannot be fully neglected in the intermediate phase containing relatively sparse oxygen ions within the VO channels.